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Abstract 



The current version of carlomat, a program for automatic computation of the lowest order 
cross sections of multiparticle reactions, is described. The program can be used as the Monte 
Carlo generator of unweighted events as well. 
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PROGRAM SUMMARY 



Title of program: carlomat 
Version: 1.0 (November 2008) 
Catalogue identifier- 
Program obtainable from: CPC Program Library or on request by e-mail from the author 
Licensing provisions: none 
Computers: all 
Operating systems: Linux 
Programming language used: FORTRAN 90/95 
Distribution format: gzipped tar archive 

Keywords: lowest order multiparticlc reactions, Standard Model, automatic calculation of 
cross sections, Monte Carlo, event generation 

Nature of physical problem 

Description of two particle scattering reactions with possibly up to 10 particles in the final 
state with a complete set of the Feynman diagrams in the lowest order of the Standard Model. 

Method of solution 

The matrix clement for a user specified process and phase space paramctrizations, which are 
necessary for the multi-channel Monte Carlo integration of the lowest order cross sections and 
event generation, are generated automatically. Both the electroweak and quantum chromody- 
namics lowest order contributions are taken into account. Particle masses are not neglected in 
the program. Matrix elements are calculated numerically with the helicity amplitude method. 
Constant widths of unstable particles are implemented by modifying mass parameters in cor- 
responding propagators. 

Restrictions on complexity of the problem 

The number of external particles is limited to 12. Only the Standard Model is implemented 

at the moment in the program. No higher order effects are taken into account, except for 
assuming the fine structure constant and the strong coupling at appropriate scale and partial 
summation of the one particle irreducible loop corrections by introducing fixed widths of 
unstable particles. 

Typical running time 

Generation of the Fortran code with carlomat on a PC with the Pentium 4 3.0 GHz processor 
for reactions with 8 particles in the final state relevant for the associated associated top 
quark pair and Higgs boson production and decay takes about 10 minutes CPU time. This 
relatively long time of the code generation is determined by a lot of write to and read from a 



disk commands which have to be introduced in order to circumvent hmitations of the Fortran 
compilers concerning possible array sizes. The compilation time of generated routines depends 
strongly on a compiler used and an optimization option chosen. Typically, for the reactions 
mentioned, it takes about one hour to compile all the routines generated. Most of the time is 
used for the compilation of the kinematical routines. The execution time of the Monte Carlo 
(MC) integration with about 2 million calls to the integrand amounts typically to a few hours 
and, if the MC summing over polarizations is employed, it is dominated by computation of 
the phase space normalization. 



LONG WRITE-UP 



1 Introduction 



Many interesting aspects of the Standard Model (SM) and models beyond it can be studied 
through investigation of reactions involving a few heavy particles at a time whose observation 
becomes possible owing to the increasing collision energy and luminosity of particle colliders, 
such as, e.g., the Large Hadron Collider (LHC), or the International Linear Collider (ILC) 
[1]. As the heavy particles are usually unstable, they almost immediately decay leading to 
reactions with several light particles in the final state which receive contributions typically 
from many thousands of Feynman diagrams already in the lowest order of SM. Although, in 
their overwhelming majority, those diagrams constitute background to the "signal diagrams" 
of production and decay of the heavy particles one wants to investigate, they must be all taken 
into account in reliable SM predictions for such reactions. This can be done practically only 
through a fully automated calculational process. 

To be more specific, let us consider, e.g. a reaction of associated production of a top quark 
pair and a Higgs boson at the ILC 

e+e" ^ tiH (1) 

that can be used to measure the top-Higgs Yukawa coupling [2]. The top and antitop decay, 
even before they hadronize, predominantly into hW^ and bW^ , respectively, the electroweak 
(EW) bosons subsequently decay into a fermion-antifermion pair each and the Higgs boson, 
if it has mass ttlh < 140 GeV, which is favoured by the direct searches at LEP and theoretical 
constrains in the framework of the SM, decays mostly into a 66-quark pair. Thus, reaction 
will actually be detected through reactions of the form 

e+e- ^ hhhhfj[hf^, (2) 

where /i,/2 = z/e, t'^, z/^, m, c and /{,/2 = , fi^ ,d, s are the decay products of the W- 
bosons coming from decays of the t- and t-quark. Except for 20 'signal' Feynman diagrams of 
the associated production of a top quark pair and a Higgs boson the representatives of which 



are depicted in Fig. [T], there are many other, typically O{10^), off resonance badcground 
diagrams which contribute to any of reactions ([2]). For example, taking into account both the 




Figure 1: Representative signal Feynman diagrams of reaction ([2]) in the unitary gauge. The 
remaining diagrams are obtained by all possible permutations of the two b and two b lines. 
The Higgs boson coupling to electrons has been neglected. 

EW and quantum chromodynamics (QCD) lowest order contributions in the unitary gauge, 
neglecting of the Yukawa couplings of the fermions lighter than c quark and r lepton, reactions 



e^e^ bbbbr^Urfi"!^^, (3) 

e^e^ bbbbud^^Ufj^, (4) 

e^e^ bbbbusdc, (5) 

e^e^ bbbbuddu (6) 



receive contributions from 21214, 26 816 and 39 342 and 185 074 Feynman diagrams, respec- 
tively. To which extent the background contributions may affect the associated production of 
the top quark pair and Higgs boson in some of the reactions listed above has been discussed 
in [3]. 

Although there exists several multipurpose Monte Carlo (MC) generators such as HELAC/ 
PHEGAS [4J, AMEGIC++/Sherpa [Sj, 'Mega/Whizard 0, MadGraph/MadEvent [7j, ALPGEN [S], 
CompHEP/CalcHEP [9j, or recently released Comix [TU], three years ago, when work on carlomat 
started, none of them was able to handle multiparticle reactions such as ([3])-([SD, not to mention 
reaction (Q. This situation has changed quite recently and now some of the publicly available 
MC generators, as e.g. O'Mega/Whizard or HELAC/ PHEGAS, are able to deliver complete SM 
predictions for reactions with 8 and more particles in the final state, including both the EW 
and QCD contributions in the leading order. However, when first physics results obtained 
with carlomat, on the off resonance background in ([T]) [3], were submitted for publication 
in spring 2008, they could only be partially checked against generators that were publicly 
available then. 

The main aim of carlomat, a new multipurpose program for automatic computation of the 
lowest order cross sections, is to provide the reliable description of multiparticle reactions. 
This is a very challenging task, not only because of complicated matrix elements, involving 
many thousands of the Feynman diagrams, but also because of the necessity of performing 
integration over a multidimensional phase space. Not only may carlomat be useful for tests 



of the complicated calculations done with other, already well established, generators. It may 
also offer the user other advantages, such as a relatively high speed of the MC computation, 
an easy control over mappings of peaks in the integrand, a possibility of calculating polarized 
cross sections, or a possibility of easy hand made modifications due to the use of traditional 
Feynman diagram approach. This write-up describes the current version of the program. 



2 Basics features of the program 

carlomat is a program written in Fortran 90/95. It generates the matrix element for a user 
specified process and phase space parametrizations which are later used for the multichannel 
MC integration of the lowest order cross sections and event generation. The program takes 
into account both the EW and QCD lowest order contributions. Particle masses are not 
neglected in the program. The number of external particles is limited to 12 and only the SM 
is currently implemented in the program. 

carlomat works according to the following scheme. The user specifies the process he/ she wants 
to have calculated. Then topologies for a given number of external particles are generated and 
checked against Feynman rules which have been coded in the program. In this process, helicity 
amplitudes, the colour matrix and phase space parametrizations are generated. Finally, they 
are copied to another directory where the numerical program can be executed. 

2.1 Generation of topologies 

Let us consider models with triple and quartic couplings. Generation of topologies starts with 
4 topologies of a 4 particle process which are depicted in Fig. |2j The 25 topologies of a 5 
particle process are obtained by attaching line No. 5 to each line, including the internal ones, 
and to each triple vertex of the graphs in Fig. [2l Particles No. 6, 7 and so on, are adeed 



recursively in the same way. The number of topologies grows dramatically with the number 
of external particles as can be seen from Table [1] However, for a process with n external 
particles, it is enough to generate topologies for n — 1 particles and then, while adding the 
n-th particle, to check whether a topology results in a Feynman diagram or not. Topologies 
can be generated and stored on a disk prior to the program execution. 




Figure 2: The topologies of a 4 particle process. 



Table 1: No. of topologies for a given No. of particles. 



No. of particles No. of topologies 



6 
7 
8 
9 
10 
11 



10525 900 
224449225 



2 485 
34 300 
559405 



220 



2.2 Feynman diagrams 

Actual initial and final state particles are assigned to lines 1, 2, 3, n in a strict order. Each 
topology is divided into two parts which are separately checked against the Feynman rules. 
While doing so, two, or three, external lines are joined by means of a triple, or quartic, vertex 
of the implemented model, respectively. In this way a new off shell particle is formed. Another 
particle is formed by adding cither the next external line, or an off shell particle that has been 
already formed. At the stage of matrix element generation, such particles are represented in 
carlomat by the derived data type called particle that is defined in the following way 

type particle 
integer : : wr,ct 
integer : : k(4) 
integer, pointer :: ip(:) 
end type particle, 



wr=l/0 if the particle has already/not yet been written to a file, 
ct=l/3/8 for a colour singlct/triplet/octet state, 

k(l) - the number of the vertex that has been used to form the particle, 

k(2) ,k(3) ,k(4) - the numbers of particles that have been joined in the vertex, with k(4)=0 

for a triple vertex, 

ip ( : ) - a list of all the off-shell particles that have been used in the definition of the current 



In the routine for computation of the matrix clement, the particles being fcrmions, vector 
bosons, or scalar bosons are represented by spinors, polarization vectors, or just complex 
numbers, respectively. 

The off shell particles and/or external particles are joined in this way until the two parts 
of a considered topology are completely covered. If they match into a propagator of the 
implemented model then the topology is accepted. Once the topology has been accepted, the 



where 



particle. 



'longer' part of it is further divided so that the Feynman diagram is made of 3 or 4 parts, joint 
to form a triple or quartic vertex of the model. This reduces the number of different particles 
that must be defined for all the Feynman diagrams. 

When the diagram is created, the corresponding particles are used to construct the helicity 
amplitude, colour factor (matrix) and phase space parametrization which are stored on the 
disk. Once all the topologies have been checked subroutines for calculating the matrix element, 
colour matrix and phase space integration are written. 



2.3 Helicity amplitudes 

Routines for computing the helicity amplitudes have been collected in a separate directory 
named carlolib. Most of the routines in carlolib had been already written for an MC 
program eett6f [11] for calculating lowest order cross sections of reactions e+e^ —>■ 6 fermions, 
relevant for tt-pair production and decay, but they have been improved and tailored to meet 
the needs of the automatic generation of amplitudes in carlomat. In order to speed up the 
computation the MC summing over helicities has been implemented in the program. To 
illustrate how it works let us consider a sum over two polarization states A = ±1 

Y: M* (A) M (A) = E ui (A) OiO.u, (A) , (7) 

A=±l A=±l 

where the polarized amplitude that depends on A has been written as 

M{\) = OaUaiX). (8) 

In Eq. ([H]), u (A) = u (p. A) is the spinor representing a particle or an antiparticle of momentum 
p and helicity A/2, O stands for the rest of amplitude M (A) that does not depend on A and a 
summation over the Dirac index a = 1, 2, 3, 4 has been assumed. Now let us define the spinor 



w [a) 



e^'^u (+1) + e-^°M (-1) , (9) 



with a being a random number uniformly distributed in the interval [0, 1]. Then the sum over 
polarizations can be replaced with the integral 

^ Mt (A) M{X) = ^ fwl {a) OlObWt (a) da. (10) 



A=±l 



Indeed, 



^ [wiia)OiO,wUa) da = 
Ztx J 



- I (e-"«i (+1) + e^-ui (-1)) OiO, {e'-u, (+1) + e'^^^ (-1)) da 



27r 





+l)OlObUk{+l)+ui{-l)OiOt,Ub{-l)= E Mt(A)M(A), 

A=±l 



27r . 27r . 2-k 

as / e^*°da = / e"^*°da = and / da = 27r. 



Analogously, a sum over three polarization states can be replaced with the integral according 
to 

M* (A) M (A) = - / e; (a) O'' t^^e, («) da, (11) 

A=0,±1 ^ 

where the longitudinal polarization component must be taken into account in the definition 
of e (a) 

e{a)=e (0) + e^^e (+1) + e^^^e (-1) (12) 

in terms of the polarization four vectors e (X) = e (p. A), A = 0, ±1, of a massive vector boson. 
Eq. ( ITT!) can be proven in exactly the same way as Eq. (fTOj) . 

Although random sampling over external helicities increases the dimension of the integral, by 
one for each helicity, it does not actually spoil its convergence, as the dependence on a in ( ITOl) 
or (fTTj) is rather smooth and the efficiency of the MC integration does not depend very much 
on the dimension of the integral. This has been explicitly checked for a number of reactions 
e^e~ — >■ 4 fermions (7) whose cross sections calculated with a program ee4f 7 [12] that applies 
explicit summing over helicities have been successfully reproduced with carlomat using the 
random sampling over external helicities. 

An explicit summing over helicities is also possible. While doing so, spinors or polarization 
vectors representing particles, both the on- and off-shell ones, are computed only once, for all 
the helicities of the external particles they are made of, and stored in arrays, which are later 
used in the sum over polarizations. 

Possible poles in the propagators of unstable particles are regularized by their constant widths 
which are introduced through the complex mass parameters 

M| = ml-imBTB, B = W,Z,H, (13) 



Mt = \/mf - imtTt, 
which replace masses in the corresponding propagators 

_ -9"' + q'qVMl 

^^^^^ " '^^TmI ' 

both in the s- and t-channels. Propagators of a photon and gluon are taken in the Feynman 
gauge. The EW mixing parameter may be defined either to be real 

sin2^^ = l-^, (14) 
mi 



with physical values of the W and Z boson masses, or complex 



sm'9w = l-j^, (15) 

z 

with and M| given by Eq. f[T51) . Using definitions ([HD, or f[T^ . in the EW coupling 
constants is usually referred to as the fixed width scheme (FWS), or the complex mass scheme 
(CMS) [13]. The colour matrix is calculated only once at the beginning of execution of the 
numerical program after having reduced its size with the use of the SU(3) algebra properties. 
The calculation is performed numerically from the very definition, using the SU(3) structure 
constants and group generators in the fundamental representation [14j . 



2.4 Phase space integration 

The integration over a multidimensional phase space is practically possible only with the use 
of Monte Carlo methods. A number of approaches has been developed in the literature to 
reduce the variance of the integrand [15]. The multichannel MC approach that is used in 
carlomat is described in this section. 

A dedicated phase space parametrization is generated for each Feynman diagram of a process 

Pi +P2 P3 + ... +p„ (16) 

with n external particles of four momenta pi, i = 1, 2, n, in the centre of mass system. The 
Lorentz invariant phase space element of (flGj) is defined in a standard way 

d'-f-'Ups = (27r)V^)L+p,-X:Jn-^, (17) 

V i=3 J i=3 W^i^i 

where rij = n — 2 is the number of particles in the final state. The set of final state particles 
is divided, in a way that depends on a topology of the diagram, into subsets consisting of two 
elements each, an element being either an external particle or a subset of particles itself and 
use is made of the identity 

/ / ^ 5^') - - = 1, n = + q,^. 

Thus, Eq. (1171) can be brought into the following form 

d^^^-Uips = (27r)^-=^'^/d/od/i...d/„_4dsids2...ds„_4, (18) 
where d/j, i = 0, 1, n — 4, is a two particle phase space element given by 



In Eq. (1191) . A is the kinematical function, g^^ and qi^ are the four momenta of each particle 
(subsystem of particles) in subset i, which are defined in their relative centre of mass system, 
+ c[j2 = 0) is the solid angle of momentum q^^ and invariants si are given by 

^. ^ f {(lii + = (^ii + EiJ^ , for i = 1, n - 4 ,^0) 
1 {Pi+ P2f = for i = 



For processes with identical particles in the initial and final states, parametrization ( 1181) must 
be slightly modified, by introducing the corresponding t invariants in a proper way. This 
has not been implemented in the program yet. Therefore, one may expect worse convergence 
of the MC integration when carlomat is used for such processes. The convergence is not 
much worse, however, if an angular cut of a few degree on the angle between the identical 
particles in the initial and final states is imposed, which has been tested for several reactions 
e^e~ 4 fermions and e^e~ 6 fermions with electrons and/or positrons in the final state. 

Invariants Si of Eq. fl20l) are randomly generated within their physical limits, s™"^ and s™^^, 
which are automatically deduced from a topology of the Feynman diagram. They can be 
generated either according to the uniform distribution 



max _ ^mm\ ^. ^ ^min^ (2I) 



where Xi denotes a random variable uniformly distributed in the interval [0, 1], or, if neces- 
sary, mappings of the Breit-Wigner shape of the propagators of unstable particles and ~ 1/s 
behaviour of the propagators of massless particles are performed. An option is included in the 
program that allows to turn on the mapping if the particle decays into 2, 3, 4, ... on-shell par- 
ticles. Different phase space parametrizations obtained in this way are important for testing 
purposes. 

Phase space parametrizations (llSp generated for each Feynman diagram 

fi{x) = &'''f-^Lipsi{x) z = l,...,iV, 

with being the number of the diagrams and x = (^Xi, x^nf-A^ uniform random arguments, 
satisfy the normalization condition 

1 

J dx3"/-V.(x) = 1. (22) 



They are used to define a new multichannel probability distribution 

TV 

f{x)=Y.aJ,{x), (23) 

i=l 

with non negative weights a^, i = 1, ...,N, satisfying the condition 



N 

J^Oi = 1 

i=l 



(24) 



which guaranties that the combined distribution f{x) is normahzed if every distribution fi{x) 
satisfies normahzation condition fl22|) . The actual MC integration is done with the random 
numbers generated according to probabihty distribution f{x). An option has been included 
in the program that allows for reducing the actual number of kinematical channels used in 
the integration. It can be turned on by choosing iopkch=l in carlocom.f. Each kinematical 
channel is then called with the same set of random arguments x and of all the weights of the 
channels with the same phase space normalization fi{x) only one weight is kept nonzero. 

The integration is performed iteratively. The first iteration starts with all the nonzero weights 
a, equal to each other, or with the weights determined in the initial scan of the integrand 
according to the following formula 

N 

ai = (Ji/^aj, (25) 
i=i 

where aj denotes the cross section obtained with the z-th kinematical channel. An option has 
been included in the program that allows to choose whether the weights should be determined 
anew according to Eq. fl25l) after every iteration, or a fraction of old weights, i.e. determined 
in the scan or previous iterations, should be preserved and transfered to the next iteration. 
This is controlled by specifying the fraction po in carlocom.f, as described in Section 3.2. 



3 Description of the program 

The two basic parts of the program which allow for a generation of the Fortran code and the 
MC computation of the cross section of a user defined process are described in the following 
two subsections. The default values of the input parameters and options used in the program 
are those specified below. 

3.1 Fortran code generation 

The part of the program that is responsible for generating routines necessary for the execution 
of the MC program is stored in a directory code_generation. The main program is stored in 
a file carlomat .f . The user should specify the process he wants to have calculated by giving 
a value to a character variable process in carlomat. f, e.g. 

process='e+ e- -> b b~mu+ vm d u~ 

He may also choose the following options in carlomat .f : 

Should the MC summing over polarizations be done, imcs=l (yes) /else (no)? 

imcs=l. Recommended for multipartide reactions. 

Take into account pure EW contributions only, iewk=l(yes)/else(no)? 

iewk=0. Recommended; takes into account both QCD and EW contributions. 



Discard the Higgs boson contributions, ihgs=l(yes)/else(no)? 

ihgs=0 . Recommended; Higgs boson contributions are taken into account. 

Then the program calls two subroutines: readpart and genmat. 

In readpart, the initial and final state particles are read from process. They are com- 
pared with the particles of the implemented model stored in a file modpart . dat and assigned 
the necessary characteristics which will be used during execution of the program. If all the 
particles match with particles in modpart . dat the possible interactions are read from a file 
vertices.dat. The numbers of vertices of different kind should not exceed the corresponding 
maximum numbers specified in a module file maxvals.f. Except for the maximum numbers 
of vertices: 

nffb=48, the number of fermion-fermion-boson vertices, 
nbbb=6, the number of triple boson vertices, 
nbbbb=8, the number of quartic boson vertices 

maxvals.f contains also values of other parameters which allow to control execution of the 
program: 

ncol=3 , the number of colors, 

mxprt=100, the maximum number of particles in a model, 
mxepl=12, the maximum number of external particles, 

mxipl=100000, the maximum number of the off-shell particles which can be created, 
mxvpl=4, the maximum number of similar vertices, such as uwy, uuZ, uuh, or uug, 
mxncnt=20 , the maximum number of continuation lines in a command, 

mxcmdl=1000 , the maximum number of lines in selected routines as, e.g., mdkk* . f , mmodkk* . f , 

mxcmdb=20, the maximum number of blocks in selected routines as, e.g., kinkk*.f , relevant 
for the speed of their compilation, 

mxmap=4, the maximum number of particles in a subsystem whose propagator should be 
mapped. 

The module contains also the name of the current version of the program. 

In genmat, topologies for the actual number of next external particles are generated, if the 
corresponding file containing prcgcncrated topologies for n = Ucxt — 1 particles could not have 
been found in a directory code_generation. This is done by consecutive calls to a subroutine 
gentop with the number of particles n = 4, 5, next- As long as n < next the topologies are 
being generated and stored in a file named topoln.dat. When gentop is called for the last 
time, with n = next, each time a new topology is created by attaching the last particle to 
a topology that was generated in the previous call and that has been now read from a file 
topoln-1 . dat, another subroutine checktop is called, which checks it against Feynman rules 
in a way described in Section 2.2. While doing so, the off-shell particles are created which, 
if the topology has been accepted, are used to generate the helicity amphtudes, the colour 



matrix and phase space parametrizations for all the Feynman diagrams corresponding to the 
topology. 

After all the topologies have been inspected, the following files necessary for a computation 
of the matrix element and phase space integration are created and stored in the working 
directory. 

• Modules: 

partcls.f, containing information about the process, as, e.g., numbers of particles, 
diagrams, etc., 

umodkk.f , a collection of complex arrays representing spinors, polarization vectors and 
scalars which are the building blocks of the diagrams, 

hmodkk . f , a collection of integer arrays which determine a sequence of elements in arrays 
of umodkk, relevant if no MC summing over polarizations is done, 

pmodkk . f , a collection of particle four momenta, 

mmodkk*.f, a collection of the amplitudes. 

• Subroutines: 

matrixel.f, the squared matrix element averaged (summed) over polarizations and 
colors of the initial (final) state particles, 

mpol2.f, the polarized squared matrix element averaged (summed) over colors of the 
initial (final) state particles, 

mtel* .f , a collection of calls to routines of carlolib which calculate building blocks of 

the Feynman diagrams, 

mdkk* .f , a collection of calls to routines of carlolib which calculate amplitudes of the 

Feynman diagrams, 

hsubkk.f , calculates integer arrays which determine a sequence of elements in arrays of 
umodkk, relevant if no MC summing over polarizations is done, 

colsqkk.f, calculates colour factors, 

kinkk*.f, a collection of mxcmdb different phase space parametrizations, 

kinds. f, calls the kinematical routines kinkk* in order to generate the final state 
particle four momenta, 

psnkk.f, calls phase space normalization routines psub*, 

psub*.f, calls the kinematical routines kinkk* in order to calculate the phase space 
normalization. 

• Parts of the code: 

partcls.dat, some information about the external particle of the process, included in 
a subroutine parf ixkk. 



mtdr*.f, adds up the amplitudes which have the same colour factors, included in a 
subroutine mpol2. 



Once the execution of carlomat ends they are copied to another directory named mc_computation, 
where the MC program can be executed. 



3.2 MC computation 

The part of the program that is responsible for the MC computation is stored in a directory 
mc_computation. The main program is stored in a file carlocom.f. In the latter, the user 
should supply the following data: 

The centre of mass energies 

aecm=(/500.d0/), 

for which the cross section should be computed together with a parameter 
ne=l 

being a size of the array aecm. If, for example, the program should run for 3 centre of mass 
energies: 500 GeV, 800 GeV and 1 TeV then aecni=(/500 .dO ,800 . dO, 1 . d3/) and ne=3. 
Moreover, the following options should be chosen in carlocom.f. 

Should the Born cross section be calculated, iborn=l(yes)/else(no) , in niter iterations 
with ncalli calls per iteration, i.e. with niter x ncalli calls? 

iborn=l, recommended, 

niter=10, 

ncalli=10000. 

The switch iborn has been introduced, since the implementation of the initial state radiation 
within the structure function approach is envisaged. Its only meaningful value in the current 
version of the program is 1. 

The fraction of old integration weights which should be transfered to the next iteration 
po=0 . dO , recommended value from to 0. 5. 

If po=0 . dO then integrations weights Oj are calculated anew according to Eq. ( l25i) and it may 
happen that, because of a low weight determined in the previous iteration, the i-th channel 
will not be selected resulting in Oj = for the next iteration. Thus the number of kinematical 
channels used in the next iterations will be effectively reduced. 

Generate the unweighted events, imc<0(no)/else(yes)? 

imc=2 . 

The probability distribution according to which the momenta of the final state particles are 
generated is limited either by the maximum of the total cross section, for imc=0, or, for 
imc=l,2, . . ., by the maxima of the cross section in each bin of distribution No. imc. See 
discussion of an option idis below. 



Should the Bom cross section be scanned with nscanO calls, iscanO=l(yes)/else(no)? 

iscanO=l , 
nscanO=10 . 

If iscanO=l then each kinematical channel is called with the same weight equal to 
where n[^°* is the total number of the channels, and the cross section is calculated with nscanO 
calls. The value of the cross section is later used to calculate the weight, according to Eq. fl25l) . 
with which that particular channel will contribute to the integral in the first iteration. 

Should the number of kinematical channels be optimized, iopkch=l(yes)/else(no)? 

iopkch=0, for processes with large number of the Feynman diagrams iopkch=l is recom- 
mended. 

Discard kinematical channels which contribute less than qw of the dominant channel in the 
initial scan. 

qw=0 . OdO . No channels are discarded; recommended. 

Choose the scheme, ischeme=l (complex mass scheme)/else(f ixed width scheme), see Eqs. 
and (|T5l) for the explanation, 

ischeme=l . 

Should distributions be calculated, idis=l(yes)/0(no)? 
idis=l . 

Variables, in which the distributions are to be calculated, the number of bins, etc., should be 
defined in a file calcdis.f. 

Should cuts be imposed, icuts=l(yes)/0(no)? 

icuts=l . 

The cuts which can be controlled from the main program carlocom currently include 

• angles: 6 (1, beam), 6 (q, beam), 6 (g, beam), 6 (a, beam), 6 (a, 1), 6 (a, q), 6 (q, q'), 6 (1, q), 

^(i,g), ^(UO 

• energies: Ex, E^, Eb, -Eg, E^ 

• invariant masses: miy, m^q', "n^bb', ""^jj', 

where 1, 1' are charged leptons, q, q' are quarks, b is a 6-quark, g is a gluon, a is a photon and 
j is either a quark or a gluon. 

The initial physical parameters: particle masses and widths, an inverse of the fine structure 
constant in the Thomson limit a, the Fermi coupling, the strong coupling constant a^, number 
of colors and the conversion constant can be specified in module inprms.f. The module 
contains one more input parameter called sqscut that is an invariant mass cut, which is used 
in mappings of the photon and gluon propagators that couple to massless fermions. Its value 



is relevant only if light fermion masses are chosen to be zero. If the width of the W boson, 
Higgs boson, or top quark is set to zero then the program will automatically replace it with a 
value calculated in the lowest order of SM. 

The helicity amplitudes are calculated with the use of the routines which have been collected 
in a directory carlolib. The directory contains also other routines as boost, kinff , kinin, 
lamsq, mapakk, mapbkk, mapskk, which are used in the phase space calculation, wwidth, a 
routine for calculating the ly-boson width in the lowest order of SM, or a random number 
generator ranlux |16] . 



4 Use of the program 

At present carlomat is distributed as a single file carlomat .tar .gz. After executing com- 
mand 

tar -xzvf carlomat . tar . gz 

a working directory named carlomat with a few subdirectories is formed. The part of the 
program responsible for code generation is stored in subdirectory code_generation. After 
having specified a desired process in carlomat . f it can be run with a command 

make test 

executed in a command line in directory code_generation. Once ended, a file test is dis- 
played on the screen that contains some details of the process under consideration, such as 
the external and off-shell particles, the number of Feynman diagrams, etc. and the generated 
files listed in Section 3.1 are moved to another working directory, named mc_computation, to 
which the user should change in order to perform the MC computation. After specifying the 
c.m.s. energy and the desired options in carlocom. f the MC program can be run by executing 
a command 

maike test. 

The output of the run is stored in a file test. Both files named test, the one in a di- 
rectory code_generation and the other in mc_computation should exactly reproduce files 
testO appended to both directories. If the MC events generation is turned on, i.e. for 
imc=0, 1,2,3, . . ., then the generated events are stored in a file mc_events.x, with x being 
the c.m.s. energy. The data for each event are stored in mc_events .x in the following form: 



# Event No. 



p( 3)= 0.47046555E+02 

p( 4)= 0.12373166E+03 

p( 5)= 0.25869511E+02 

p( 6)= 0.17716225E+03 

p( 7)= 0.65263801E+02 

p( 8)= 0.60926224E+02 



-0.10891034E+02 
0.59317748E+02 
0.21415118E+02 
-0.40997709E+02 
-0.52555950E+02 
0.23711829E+02 



0.42906506E+02 
0.19678931E+02 
-0.37079632E+01 
0.81425324E+01 
-0.20628053E+02 
-0.46391953E+02 



0.15190637E+02 
0.10668003E+03 
-0.14030830E+02 
-0.17216082E+03 
0.32737122E+02 
0.31583866E+02 



where p(3) , . . . ,p(8) are four momenta of the final state particles in the centre of mass 
system in exactly the same order as they appear in process. The first component is the 
energy and the other 3 components are the x, y and z components of the four momentum. 
The number of accepted events substantially increases when the program is run the second 
time with exactly the same choice of options and parameters. In each run, the maximum 
values of the cross section are stored in a file named csmax_* in directory mc_computation 
and they are automatically read from the disc in the next run of the program. The efficiency of 
events acceptance is not high, but it can be improved by a proper selection of the differential 
cross section and an appropriate choice of the number of bins. 



5 Outlook 



The following improvements of the program are envisaged in the near future. 

• Dedicated treatment of soft and collinear external photons, as well as t-channel pho- 
ton/gluon exchange. 

• Interfaces to parton density functions, or the initial state radiation within the structure 
function approach. 

• More efficient calculation of colour coefficients. 

Moreover, interfaces to parton shower and hadronization programs should be worked on. Ex- 
tensions of SM can be implemented and the corresponding lowest order cross sections can be 
calculated in a fully automatic way. Leading SM radiative corrections can be implemented, if 
the necessary subroutines are provided. 
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